In [142]:
%pylab inline
from scipy.optimize import basinhopping
from scipy.interpolate import UnivariateSpline
import pandas as pd
import csv
import shapely.affinity as af
import shapely.geometry as sh
from equivalent_ellipse import *
from scaled_figures import *
from descartes.patch import PolygonPatch
In [143]:
pylab.rcParams['savefig.dpi'] = 110
In [144]:
x_list = list()
y_list = list()
with open('../data/custom_cutout_x.csv', 'r') as x_csvfile:
with open('../data/custom_cutout_y.csv', 'r') as y_csvfile:
x_reader = csv.reader(x_csvfile, delimiter=',', lineterminator='\n')
y_reader = csv.reader(y_csvfile, delimiter=',', lineterminator='\n')
for row in x_reader:
x_list += [row]
for row in y_reader:
y_list += [row]
num_cutouts = len(x_list)
x_array = [0,]*num_cutouts
y_array = [0,]*num_cutouts
for i in range(num_cutouts):
x_array[i] = array(x_list[i], dtype='float')
y_array[i] = array(y_list[i], dtype='float')
cutout = [0,]*num_cutouts
for i in range(num_cutouts):
cutout[i] = shapely_cutout(x_array[i],y_array[i])
In [145]:
cutout[0]
Out[145]:
In [146]:
# Use pandas to load the csv file
cutoutData = pd.read_csv("cutoutdata.csv")
# Define width, length, and factor as the csv columns
# which had that corresponding header
width = cutoutData['width'].values
length = cutoutData['length'].values
factor = cutoutData['factor'].values
cutoutRef = cutoutData['Reference'].values
# Swap width and length values if order was mistaken
swap = width > length
store = width[swap]
width[swap] = length[swap]
length[swap] = store
In [147]:
cutoutRef = cutoutRef[~isnan(cutoutRef)].astype(int)
cutoutRef
Out[147]:
In [148]:
# factor
In [149]:
# with open('cutout_factors','r') as f:
# loaded_factor = eval(f.read())
# dictArray = array(list(loaded_factor.items()),float)
# cutoutRef = dictArray[:,0].astype(int)
# factor = dictArray[:,1]
# # width = width[cutoutRef]
# # length = length[cutoutRef]
# factor
In [150]:
# with open('custom_cutout_factors','r') as f:
# loaded_custom_factor = eval(f.read())
# customDimensionsDict = {'5x13':[5,13],
# '5x10':[5,10],
# '5x8':[5,8],
# '4x13':[4,13],
# '4x10':[4,10],
# '4x8':[4,8],
# '4x6.5':[4,6.5],
# '3x13':[3,13],
# '3x9':[3,9],
# '3x6.5':[3,6.5],
# '3x5':[3,5]}
# # customDimensionsDict = {'5x8':[5,8],
# # '4x8':[4,8],
# # '4x6.5':[4,6.5],
# # '3x9':[3,9],
# # '3x6.5':[3,6.5],
# # '3x5':[3,5]}
# customDataLabelList = list(customDimensionsDict.keys())
# custom_factors = [loaded_custom_factor[x] for x in customDataLabelList]
# custom_widths = [customDimensionsDict[x][0] for x in customDataLabelList]
# custom_lengths = [customDimensionsDict[x][1] for x in customDataLabelList]
# ellipse_cutouts = [0,]*len(custom_widths)
# for i in range(len(custom_widths)):
# ellipse_cutouts[i] = create_ellipse([0,0,custom_widths[i],custom_lengths[i],-45])
# cutoutRef = append(cutoutRef,arange(len(cutout),len(cutout)+len(ellipse_cutouts)))
# cutout.extend(ellipse_cutouts)
# factor = append(factor,custom_factors)
# width = append(width,custom_widths)
# length = append(length,custom_lengths)
In [151]:
# testCutout = cutout[70]
# testCutout
In [152]:
with open('circle_cutout_factors','r') as f:
loaded_circle_factor = eval(f.read())
circleDictArray = array(list(loaded_circle_factor.items()),float)
circleDiameter = circleDictArray[:,0]
circleFactor = circleDictArray[:,1]
circle_cutouts = [0,]*len(circleDiameter)
for i in range(len(circleDiameter)):
circle_cutouts[i] = create_ellipse([0,0,circleDiameter[i],circleDiameter[i],0])
In [153]:
scatter(width,length)
Out[153]:
In [154]:
scatter(width, width/length)
Out[154]:
In [155]:
scatter(width, 1/factor)
circleFit = UnivariateSpline(circleDiameter,1/circleFactor)
diameter_i = linspace(circleDiameter.min(),circleDiameter.max())
plot(diameter_i,circleFit(diameter_i),'r--')
widthFit = UnivariateSpline(width,1/factor)
width_i = linspace(width.min(),width.max())
plot(width_i,widthFit(width_i),'k--')
Out[155]:
In [156]:
pylab.rcParams['savefig.dpi'] = 90
In [157]:
ratio = width/length
artificial = isnan(cutoutData['Reference'].values)
scatter(width[~artificial],ratio[~artificial])
circleRef = ratio == 1
scatter(width[circleRef],ratio[circleRef],c='red')
scatter(width[artificial & ~circleRef],ratio[artificial & ~circleRef],c='green')
title("Distrbution of selected cutouts")
xlabel("Width (cm)")
ylabel("Ratio")
legend(['Clinical','Circles','Ellipses'],loc=4)
Out[157]:
In [158]:
scatter(width[~artificial],1/factor[~artificial])
circleRef = ratio == 1
scatter(width[circleRef],1/factor[circleRef],c='red')
scatter(width[artificial & ~circleRef],1/factor[artificial & ~circleRef],c='green')
# widthFit = UnivariateSpline(width,1/factor)
# width_i = linspace(width.min(),width.max())
plot(width_i,widthFit(width_i),'r--')
title("The effect of width on cutout factor")
ylabel("1 / cutout factor")
xlabel("Width cm")
Out[158]:
In [159]:
# scatter(width, 1/factor)
pylab.rcParams['savefig.dpi'] = 90
# widthFit = UnivariateSpline(width,1/factor)
width_i = linspace(width.min(),width.max(),1000)
derivFunc = widthFit.derivative()
deriv = widthFit.derivative()(width_i)
deriv[deriv<0] = 0
plot(width_i,deriv,'b')
plot(width_i,zeros(len(width_i)),'k--')
# scatter(width_i,deriv,c=deriv,edgecolors='none',s=2)
title("Contribution to cutout factor\nfrom ring of given diameter")
ylabel("Derivative of (1 / cutout factor wrt width)")
xlabel("Diameter | Width (cm)")
xlim([3,9])
ylim([-0.002,0.025])
Out[159]:
In [159]:
In [160]:
minVal = 1.5
maxVal = 15
numZones = 1000
In [161]:
# Create Zones
zoneWidth = linspace(minVal,maxVal,numZones)
zoneWeight = 70*widthFit.derivative()(zoneWidth)
zoneWeight[zoneWeight < 0] = 0
circle = [0,]*(numZones)
circle[0] = sh.Point(0,0).buffer(minVal/2)
zone = [0,]*numZones
zone[0] = circle[0]
for i in range(1,numZones-1):
circle[i] = sh.Point(0,0).buffer(zoneWidth[i]/2)
zone[i] = circle[i].difference(circle[i-1])
zone[-1] = sh.Point(0,0).buffer(20).difference(circle[-2])
weightedArea = 0.
for i in range(numZones):
weightedArea += cutout[0].intersection(zone[i]).area * zoneWeight[i]
weightedArea
Out[161]:
In [162]:
cutout[0].area
Out[162]:
In [163]:
plot(zoneWidth,zoneWeight)
Out[163]:
In [164]:
# def weightedAreaOptimise(x,inputShape):
# a = x[0]
# b = x[1]
# movedShape =af.translate(inputShape, xoff=-a, yoff=-b)
# weightedArea = 0.
# for i in range(numZones):
# weightedArea += movedShape.intersection(zone[i]).area * zoneWeight[i]
# return -weightedArea
In [165]:
# class MyTakeStep(object):
# def __init__(self, stepsize=0.5):
# self.stepsize = stepsize
# def __call__(self, x):
# s = self.stepsize
# x[0] += np.random.uniform(-3*s, 3*s)
# x[1] += np.random.uniform(-3*s, 3*s)
# return x
# def print_fun(x, f, accepted):
# global successCount_basinhopping
# global minVals_basinhopping
# global nSuccess_basinhopping
# a = x[0]
# b = x[1]
# f = np.around(f,4)
# if accepted:
# if np.all(np.isnan(minVals_basinhopping)):
# minVals_basinhopping[0] = f
# successCount_basinhopping = 1
# print(("first local minima of %.4f at:\n"+
# "[a,b] = [%.4f, %.4f]")
# % (f,a,b))
# elif (np.nanmin(minVals_basinhopping) == f):
# minVals_basinhopping[successCount_basinhopping] = f
# successCount_basinhopping += 1
# print(("agreeing local minima of %.4f at:\n"+
# "[a,b] = [%.4f, %.4f]")
# % (f,a,b))
# elif (np.nanmin(minVals_basinhopping) > f):
# minVals_basinhopping[0] = f
# successCount_basinhopping = 1
# print(("new local minima of %.4f at:\n"+
# "[a,b] = [%.4f, %.4f]")
# % (f,a,b))
# else:
# print(("rejected local minima of %.4f at:\n"+
# "[a,b] = [%.4f, %.4f]")
# % (f,a,b))
# else:
# print(("failed to find local minima at:\n"+
# "[a,b] = [%.4f, %.4f]")
# % (a,b))
# if successCount_basinhopping >= nSuccess_basinhopping:
# return True
# def optimise_middle(cutout,userNSuccess):
# mytakestep = MyTakeStep()
# global nSuccess_basinhopping
# nSuccess_basinhopping = userNSuccess
# global minVals_basinhopping
# minVals_basinhopping = np.empty(nSuccess_basinhopping)
# minVals_basinhopping[:] = np.nan
# global successCount_basinhopping
# successCount_basinhopping = 0
# minimizer_kwargs = {"args": (cutout,), "method": 'L-BFGS-B',
# "bounds": ((-5,5),(-5,5)) }
# x0 = np.array([0,0])
# output = basinhopping(weightedAreaOptimise,x0,
# niter=1000,minimizer_kwargs=minimizer_kwargs,
# take_step=mytakestep, callback=print_fun)
# return output
In [166]:
# minVal = 2
# maxVal = 8
# numZones = 100
# zoneWidth = linspace(minVal,maxVal,numZones)
# zoneWeight = 70*widthFit.derivative()(zoneWidth)
# zoneWeight[zoneWeight < 0] = 0
# circle = [0,]*(numZones)
# circle[0] = sh.Point(0,0).buffer(minVal/2)
# zone = [0,]*numZones
# zone[0] = circle[0]
# for i in range(1,numZones-1):
# circle[i] = sh.Point(0,0).buffer(zoneWidth[i]/2)
# zone[i] = circle[i].difference(circle[i-1])
# zone[-1] = sh.Point(0,0).buffer(20).difference(circle[-2])
# def weightedAreaCalc(inputShape):
# weightedArea = 0.
# for i in range(numZones):
# weightedArea += inputShape.intersection(zone[i]).area * zoneWeight[i]
# return weightedArea
# def weightedAreaDifference(cutout,ellipse,a,b):
# cutout = af.translate(cutout, xoff=-a, yoff=-b)
# ellipse = af.translate(ellipse, xoff=-a, yoff=-b)
# # inside_factor = ellipse.difference(cutout).area / ellipse.area
# # inside_factor = exp(-(inside_factor-0.2))**40
# # print(inside_factor)
# return (weightedAreaCalc(cutout) +
# weightedAreaCalc(ellipse) -
# weightedAreaCalc(ellipse.difference(cutout)) -
# weightedAreaCalc(cutout.difference(ellipse)))
In [167]:
# def circle_in_cutout(minimiserInput,cutout):
# x = zeros(5)
# x[0] = minimiserInput[0]
# x[1] = minimiserInput[1]
# x[2] = minimiserInput[2]
# x[3] = minimiserInput[2]
# x[4] = 0
# circle = create_ellipse(x)
# toMinimise = -100*weightedAreaDifference(cutout,circle,x[0],x[1])
# return toMinimise
# class MyTakeStep(object):
# def __init__(self, stepsize=0.5):
# self.stepsize = stepsize
# def __call__(self, x):
# s = self.stepsize
# x[0] += np.random.uniform(-3*s, 3*s)
# x[1] += np.random.uniform(-3*s, 3*s)
# x[2] += np.random.uniform(-6*s, 6*s)
# return x
# def print_fun(x, f, accepted):
# global successCount_basinhopping
# global minVals_basinhopping
# global nSuccess_basinhopping
# a = x[0]
# b = x[1]
# f = np.around(f,4)
# if accepted:
# if np.all(np.isnan(minVals_basinhopping)):
# minVals_basinhopping[0] = f
# successCount_basinhopping = 1
# print(("first local minima of %.4f at:\n"+
# "[a,b] = [%.4f, %.4f]")
# % (f,a,b))
# elif (np.nanmin(minVals_basinhopping) == f):
# minVals_basinhopping[successCount_basinhopping] = f
# successCount_basinhopping += 1
# print(("agreeing local minima of %.4f at:\n"+
# "[a,b] = [%.4f, %.4f]")
# % (f,a,b))
# elif (np.nanmin(minVals_basinhopping) > f):
# minVals_basinhopping[0] = f
# successCount_basinhopping = 1
# print(("new local minima of %.4f at:\n"+
# "[a,b] = [%.4f, %.4f]")
# % (f,a,b))
# else:
# print(("rejected local minima of %.4f at:\n"+
# "[a,b] = [%.4f, %.4f]")
# % (f,a,b))
# else:
# print(("failed to find local minima at:\n"+
# "[a,b] = [%.4f, %.4f]")
# % (a,b))
# if successCount_basinhopping >= nSuccess_basinhopping:
# return True
# def optimise_middle(cutout,userNSuccess):
# mytakestep = MyTakeStep()
# global nSuccess_basinhopping
# nSuccess_basinhopping = userNSuccess
# global minVals_basinhopping
# minVals_basinhopping = np.empty(nSuccess_basinhopping)
# minVals_basinhopping[:] = np.nan
# global successCount_basinhopping
# successCount_basinhopping = 0
# minimizer_kwargs = {"args": (cutout,), "method": 'L-BFGS-B',
# "bounds": ((-5,5),(-5,5),
# (minVal,maxVal)) }
# x0 = np.array([0,0,3])
# output = basinhopping(circle_in_cutout,x0,
# niter=1000,minimizer_kwargs=minimizer_kwargs,
# take_step=mytakestep, callback=print_fun)
# # output = basinhopping(circle_in_cutout,x0,
# # niter=1000,minimizer_kwargs=minimizer_kwargs)
# return output
In [168]:
# zoneWidth = linspace(minVal,maxVal,numZones)
# zoneWeight = 70*widthFit.derivative()(zoneWidth)
# zoneWeight[zoneWeight < 0] = 0
# circle = [0,]*(numZones)
# circle[0] = sh.Point(0,0).buffer(minVal/2)
# zone = [0,]*numZones
# zone[0] = circle[0]
# for i in range(1,numZones-1):
# circle[i] = sh.Point(0,0).buffer(zoneWidth[i]/2)
# zone[i] = circle[i].difference(circle[i-1])
# zone[-1] = sh.Point(0,0).buffer(20).difference(circle[-2])
def weightedAreaCalc(inputShape):
weightedArea = 0.
for i in range(numZones):
weightedArea += inputShape.intersection(zone[i]).area * zoneWeight[i]
return weightedArea
def weightedAreaDifference(cutout,ellipse,a,b):
cutout = af.translate(cutout, xoff=-a, yoff=-b)
ellipse = af.translate(ellipse, xoff=-a, yoff=-b)
return (weightedAreaCalc(ellipse.difference(cutout)) +
weightedAreaCalc(cutout.difference(ellipse)))
# def weightedAreaDifference(cutout,ellipse,a,b):
# cutout = af.translate(cutout, xoff=-a, yoff=-b)
# ellipse = af.translate(ellipse, xoff=-a, yoff=-b)
# return ellipse.difference(cutout).area + cutout.difference(ellipse).area
def ellipse_in_cutout(minimiserInput,cutout,weight):
ellipse = create_ellipse(minimiserInput)
a = minimiserInput[0]
b = minimiserInput[1]
toMinimise = 100*weightedAreaDifference(cutout,ellipse,a,b)
return toMinimise
class MyTakeStep(object):
def __init__(self, stepsize=0.5):
self.stepsize = stepsize
def __call__(self, x):
s = self.stepsize
x[0] += np.random.uniform(-3*s, 3*s)
x[1] += np.random.uniform(-3*s, 3*s)
x[2] += np.random.uniform(-6*s, 6*s)
x[3] += np.random.uniform(-8*s, 8*s)
x[4] += np.random.uniform(-90*s, 90*s)
return x
def print_fun(x, f, accepted):
global successCount_basinhopping
global minVals_basinhopping
global nSuccess_basinhopping
a = x[0]
b = x[1]
w = abs(x[2])
l = abs(x[3])
if (l > w):
width = w
length = l
theta = np.mod(x[4],180)
else:
width = l
length = w
theta = np.mod(x[4]+90,180)
f = np.around(f,4)
if accepted:
if np.all(np.isnan(minVals_basinhopping)):
minVals_basinhopping[0] = f
successCount_basinhopping = 1
print(("first local minima of %.4f at:\n"+
"[a,b] = [%.4f, %.4f],"+
" [width, length] = [%.4f, %.4f],"+
" theta = %.4f\n")
% (f,a,b,width,length,theta))
elif (np.nanmin(minVals_basinhopping) == f):
minVals_basinhopping[successCount_basinhopping] = f
successCount_basinhopping += 1
print(("agreeing local minima of %.4f at:\n"+
"[a,b] = [%.4f, %.4f],"+
" [width, length] = [%.4f, %.4f],"+
" theta = %.4f\n")
% (f,a,b,width,length,theta))
elif (np.nanmin(minVals_basinhopping) > f):
minVals_basinhopping[0] = f
successCount_basinhopping = 1
print(("new local minima of %.4f at:\n"+
"[a,b] = [%.4f, %.4f],"+
" [width, length] = [%.4f, %.4f],"+
" theta = %.4f\n")
% (f,a,b,width,length,theta))
else:
print(("rejected local minima of %.4f at:\n"+
"[a,b] = [%.4f, %.4f],"+
" [width, length] = [%.4f, %.4f],"+
" theta = %.4f\n")
% (f,a,b,width,length,theta))
else:
print(("failed to find local minima at:\n"+
"[a,b] = [%.4f, %.4f],"+
" [width, length] = [%.4f, %.4f],"+
" theta = %.4f\n")
% (a,b,width,length,theta))
if successCount_basinhopping >= nSuccess_basinhopping:
return True
def optimise_middle(cutout,weight,userNSuccess):
mytakestep = MyTakeStep()
global nSuccess_basinhopping
nSuccess_basinhopping = userNSuccess
global minVals_basinhopping
minVals_basinhopping = np.empty(nSuccess_basinhopping)
minVals_basinhopping[:] = np.nan
global successCount_basinhopping
successCount_basinhopping = 0
minimizer_kwargs = {"args": (cutout,weight), "method": 'L-BFGS-B',
"bounds": ((-5,5),(-5,5),
(minVal,maxVal),
(minVal,maxVal),
(None, None)) }
x0 = np.array([0,0,3,5,0])
output = basinhopping(ellipse_in_cutout,x0,
niter=1000,minimizer_kwargs=minimizer_kwargs,
take_step=mytakestep, callback=print_fun)
return output
In [169]:
ellipse_results = [0,]*len(cutoutRef)
# for j in range(len(cutoutRef)):
for j in [0]:
i = j
print("Cutout %.f" % (i))
ellipse_results[j] = optimise_middle(cutout[i],1,2)
print(" ")
In [170]:
ellipse_def = zeros([len(cutoutRef),5])
# for j in range(len(cutoutRef)):
for j in [0]:
ellipse_def[j,:] = array(ellipse_results[j].x)
In [171]:
middle_points = zeros([len(cutoutRef),2])
# for j in range(len(cutoutRef)):
for j in [0]:
i = j
middle_points[j,:] = ellipse_def[j,0:2]
In [172]:
# j = find(cutoutRef == 14)[0]
# i = cutoutRef[j]
i = 0; j = 0
ellipse = create_ellipse(ellipse_def[j,:])
scaled_fig_start(12,12)
# plot(centred_cutout[j].exterior.xy[0],centred_cutout[j].exterior.xy[1])
plot(cutout[i].exterior.xy[0],cutout[i].exterior.xy[1],'b')
plot(ellipse.exterior.xy[0],ellipse.exterior.xy[1],'g',linewidth=4,alpha=0.5)
scatter(middle_points[j,0],middle_points[j,1],c='g')
text(-5,5,('Width: %.2f cm' % (ellipse_def[j,2])),fontsize=11)
text(-5,4.5,('Length: %.2f cm' % (ellipse_def[j,3])),fontsize=11)
scaled_fig_end(12,12)
In [172]:
In [172]:
In [172]:
In [173]:
# x = zeros(5)
# x[0] = output.x[0]
# x[1] = output.x[1]
# x[2] = output.x[2]
# x[3] = output.x[2]
# x[4] = 0
# circle = create_ellipse(x)
In [174]:
# circle.intersection(cutout[0]).area
In [175]:
# j = 0
# i = 0
# a,b = middle_points[j,:]
# scaled_fig_start(10,10)
# plot(cutout[i].exterior.xy[0],cutout[i].exterior.xy[1])
# plot(circle.exterior.xy[0],circle.exterior.xy[1])
# plot(a,b,'go')
# scaled_fig_end(10,10)
In [176]:
testcutout = af.translate(cutout[i], xoff=-a, yoff=-b)
testcutout
Out[176]:
In [177]:
minVal = 1.5
maxVal = 15
numZones = 1000
# Create Zones
zoneWidth = linspace(minVal,maxVal,numZones)
zoneWeight = 70*widthFit.derivative()(zoneWidth)
zoneWeight[zoneWeight < 0] = 0
circle = [0,]*(numZones)
circle[0] = sh.Point(0,0).buffer(minVal/2)
zone = [0,]*numZones
zone[0] = circle[0]
for i in range(1,numZones-1):
circle[i] = sh.Point(0,0).buffer(zoneWidth[i]/2)
zone[i] = circle[i].difference(circle[i-1])
zone[-1] = sh.Point(0,0).buffer(20).difference(circle[-2])
In [178]:
ratioAreaInZones = zeros(numZones)
for k in range(numZones):
ratioAreaInZones[k] = testcutout.intersection(zone[k]).area / zone[k].area
plot(zoneWidth,ratioAreaInZones)
Out[178]:
In [179]:
dotAngle = 90-90*ratioAreaInZones
plot(zoneWidth,dotAngle)
Out[179]:
In [180]:
x_store = array([])
y_store = array([])
for k in range(len(dotAngle)):
if (dotAngle[k] > 0) & (dotAngle[k] < 90):
hypot = zoneWidth[k]/2
y_val = hypot * sin(pi/180*dotAngle[k])
x_val = hypot * cos(pi/180*dotAngle[k])
x_store = append(x_store,x_val)
y_store = append(y_store,y_val)
x_mirror = array([x_val,x_val,-x_val,-x_val])
y_mirror = array([y_val,-y_val,y_val,-y_val])
scatter(x_mirror,y_mirror)
axis('equal')
Out[180]:
In [181]:
ordered_x = concatenate((x_store,-x_store[::-1],-x_store,x_store[::-1]))
ordered_y = concatenate((y_store,y_store[::-1],-y_store,-y_store[::-1]))
In [182]:
straightCutout = shapely_cutout(ordered_x,ordered_y)
# straightCutout = af.rotate(straightCutout, -45)
straightCutout
Out[182]:
In [199]:
def visualRotation(x,original,new):
rotatedShape = af.rotate(new, x)
shapeIntersection = rotatedShape.intersection(original)
# weightedOverlap = 0.
# for i in range(numZones):
# weightedOverlap += shapeIntersection.intersection(zone[i]).area * zoneWeight[i]
# print(-shapeIntersection.area)
return -shapeIntersection.area
minimizer_kwargs = {"args": (testcutout,straightCutout,), "method": 'Powell'}
x0 = np.array([-90])
# output = basinhopping(visualRotation,x0,
# niter=3,minimizer_kwargs=minimizer_kwargs)
straightCutout = af.rotate(straightCutout, -45)
output
Out[199]:
In [200]:
scaled_fig_start(10,10)
plot(testcutout.exterior.xy[0],testcutout.exterior.xy[1],linewidth=2)
plot(straightCutout.exterior.xy[0],straightCutout.exterior.xy[1],'k',linewidth=2,alpha=1)
scatter(0,0,c='green',s=30)
scaled_fig_end(10,10)
In [201]:
scaled_fig_start(10,10)
# plot(straightCutout.exterior.xy[0],straightCutout.exterior.xy[1],'k',linewidth=2,alpha=0.5)
scatter(0,0,c='green',s=40)
for k in range(0,numZones,100):
plot(zone[k].exterior.xy[0],zone[k].exterior.xy[1],'k',linewidth=2,alpha=0.5)
plot(testcutout.exterior.xy[0],testcutout.exterior.xy[1],linewidth=2)
overlap = zone[2]
scaled_fig_end(10,10)
In [202]:
# scatter(width, 1/factor)
# widthFit = UnivariateSpline(width,1/factor)
width_i = linspace(width.min(),15,1000)
derivFunc = widthFit.derivative()
deriv = widthFit.derivative()(width_i)
# plot(width_i,deriv,'b')
plot(width_i,zeros(len(width_i)),'k--')
scatter(width_i,deriv,c=deriv,edgecolors='none',s=2)
title("Contribution to cutout factor\nfrom ring of given diameter")
ylabel("Derivative of (1 / cutout factor wrt width)")
xlabel("Diameter | Width (cm)")
# xlim([3,15])
# ylim([-0.005,0.025])
Out[202]:
In [203]:
numRands = 100000
x0 = 0
y0 = 0
dist = width.min()/2 + rand(numRands) * width.ptp()/2
angle = 2*pi*rand(numRands)
xVals = x0 + dist*cos(angle)
yVals = y0 + dist*sin(angle)
scaled_fig_start(10,10)
scatter(xVals,yVals,c=derivFunc(dist),edgecolors='none',s=5,alpha=0.1)
plot(testcutout.exterior.xy[0],testcutout.exterior.xy[1],'k',linewidth=2)
scaled_fig_end(10,10)
In [203]:
In [204]:
scaled_fig_start(10,10)
# plot(straightCutout.exterior.xy[0],straightCutout.exterior.xy[1],'k',linewidth=2,alpha=0.5)
scatter(0,0,c='green',s=40)
for k in range(0,numZones,100):
plot(zone[k].exterior.xy[0],zone[k].exterior.xy[1],'k',linewidth=1,alpha=0.5)
plot(testcutout.exterior.xy[0],testcutout.exterior.xy[1],linewidth=2)
overlap = zone[2]
scaled_fig_end(10,10)
In [205]:
numZonesTest = 10
zoneWidthTest = linspace(minVal,maxVal,numZonesTest)
# zoneWeight = 70*widthFit.derivative()(zoneWidth)
circleTest = [0,]*(10)
circleTest[0] = sh.Point(0,0).buffer(minVal/2)
zoneTest = [0,]*numZonesTest
zoneTest[0] = circleTest[0]
for i in range(1,numZonesTest-1):
circleTest[i] = sh.Point(0,0).buffer(zoneWidthTest[i]/2)
zoneTest[i] = circleTest[i].difference(circleTest[i-1])
zoneTest[-1] = sh.Point(0,0).buffer(20).difference(circleTest[-2])
In [206]:
pylab.rcParams['savefig.dpi'] = 70
BLUE = '#6699cc'
YELLOW = '#ffcc33'
GREEN = '#339933'
GRAY = '#999999'
RED = '#ff3333'
In [207]:
scaled_fig_start(10,10)
ax = gca()
# plot(straightCutout.exterior.xy[0],straightCutout.exterior.xy[1],'k',linewidth=2,alpha=0.5)
scatter(0,0,c='green',s=40)
for k in range(0,numZonesTest):
plot(zoneTest[k].exterior.xy[0],zoneTest[k].exterior.xy[1],'k',linewidth=1,alpha=0.5)
intersect = testcutout.intersection(zoneTest[2])
if intersect.geom_type == 'Polygon':
intersectPatch = PolygonPatch(intersect, fc=GREEN, ec=GREEN, alpha=0.5, zorder=1)
ax.add_patch(intersectPatch)
elif intersect.geom_type == 'MultiPolygon':
for p in intersect:
intersectPatch = PolygonPatch(p, fc=GREEN, ec=GREEN, alpha=0.5, zorder=1)
ax.add_patch(intersectPatch)
# overlap = PolygonPatch(straightCutout.intersection(zoneTest[3]), fc=BLUE, ec=BLUE, alpha=0.5, zorder=1)
# ax.add_patch(overlap)
print("Area = %0.1f cm^2"% (intersect.area,))
intersect = testcutout.intersection(zoneTest[3])
if intersect.geom_type == 'Polygon':
intersectPatch = PolygonPatch(intersect, fc=BLUE, ec=BLUE, alpha=0.5, zorder=1)
ax.add_patch(intersectPatch)
elif intersect.geom_type == 'MultiPolygon':
for p in intersect:
intersectPatch = PolygonPatch(p, fc=BLUE, ec=BLUE, alpha=0.5, zorder=1)
ax.add_patch(intersectPatch)
# overlap = PolygonPatch(straightCutout.intersection(zoneTest[3]), fc=BLUE, ec=BLUE, alpha=0.5, zorder=1)
# ax.add_patch(overlap)
print("Area = %0.1f cm^2"% (intersect.area,))
intersect = testcutout.intersection(zoneTest[5])
if intersect.geom_type == 'Polygon':
intersectPatch = PolygonPatch(intersect, fc=RED, ec=RED, alpha=0.5, zorder=1)
ax.add_patch(intersectPatch)
elif intersect.geom_type == 'MultiPolygon':
for p in intersect:
intersectPatch = PolygonPatch(p, fc=RED, ec=RED, alpha=0.5, zorder=1)
ax.add_patch(intersectPatch)
# overlap = PolygonPatch(straightCutout.intersection(zoneTest[3]), fc=BLUE, ec=BLUE, alpha=0.5, zorder=1)
# ax.add_patch(overlap)
print("Area = %0.1f cm^2"% (intersect.area,))
plot(testcutout.exterior.xy[0],testcutout.exterior.xy[1],linewidth=2)
scaled_fig_end(10,10)
In [208]:
scaled_fig_start(10,10)
ax = gca()
# plot(straightCutout.exterior.xy[0],straightCutout.exterior.xy[1],'k',linewidth=2,alpha=0.5)
scatter(0,0,c='green',s=40)
for k in range(0,numZonesTest):
plot(zoneTest[k].exterior.xy[0],zoneTest[k].exterior.xy[1],'k',linewidth=1,alpha=0.5)
intersect = straightCutout.intersection(zoneTest[2])
if intersect.geom_type == 'Polygon':
intersectPatch = PolygonPatch(intersect, fc=GREEN, ec=GREEN, alpha=0.5, zorder=1)
ax.add_patch(intersectPatch)
elif intersect.geom_type == 'MultiPolygon':
for p in intersect:
intersectPatch = PolygonPatch(p, fc=GREEN, ec=GREEN, alpha=0.5, zorder=1)
ax.add_patch(intersectPatch)
# overlap = PolygonPatch(straightCutout.intersection(zoneTest[3]), fc=BLUE, ec=BLUE, alpha=0.5, zorder=1)
# ax.add_patch(overlap)
print("Area = %0.1f cm^2"% (intersect.area,))
intersect = straightCutout.intersection(zoneTest[3])
if intersect.geom_type == 'Polygon':
intersectPatch = PolygonPatch(intersect, fc=BLUE, ec=BLUE, alpha=0.5, zorder=1)
ax.add_patch(intersectPatch)
elif intersect.geom_type == 'MultiPolygon':
for p in intersect:
intersectPatch = PolygonPatch(p, fc=BLUE, ec=BLUE, alpha=0.5, zorder=1)
ax.add_patch(intersectPatch)
# overlap = PolygonPatch(straightCutout.intersection(zoneTest[3]), fc=BLUE, ec=BLUE, alpha=0.5, zorder=1)
# ax.add_patch(overlap)
print("Area = %0.1f cm^2"% (intersect.area,))
intersect = straightCutout.intersection(zoneTest[5])
if intersect.geom_type == 'Polygon':
intersectPatch = PolygonPatch(intersect, fc=RED, ec=RED, alpha=0.5, zorder=1)
ax.add_patch(intersectPatch)
elif intersect.geom_type == 'MultiPolygon':
for p in intersect:
intersectPatch = PolygonPatch(p, fc=RED, ec=RED, alpha=0.5, zorder=1)
ax.add_patch(intersectPatch)
# overlap = PolygonPatch(straightCutout.intersection(zoneTest[3]), fc=BLUE, ec=BLUE, alpha=0.5, zorder=1)
# ax.add_patch(overlap)
print("Area = %0.1f cm^2"% (intersect.area,))
plot(straightCutout.exterior.xy[0],straightCutout.exterior.xy[1],linewidth=2)
scaled_fig_end(10,10)
In [209]:
centred_cutout = [0.]*len(cutoutRef)
straightened_cutout = [0.]*len(cutoutRef)
visual_angle = zeros(len(cutoutRef))
# for j in range(len(cutoutRef)):
for j in [0]:
i = j
print("Cutout %d" %i)
a,b = middle_points[j,:]
centred_cutout[j] = af.translate(cutout[i], xoff=-a, yoff=-b)
ratioAreaInZones = zeros(numZones)
for k in range(numZones):
ratioAreaInZones[k] = centred_cutout[j].intersection(zone[k]).area / zone[k].area
dotAngle = 90-90*ratioAreaInZones
x_store = array([])
y_store = array([])
for k in range(len(dotAngle)):
if (dotAngle[k] > 0) & (dotAngle[k] < 90):
hypot = zoneWidth[k]/2
y_val = hypot * sin(pi/180*dotAngle[k])
x_val = hypot * cos(pi/180*dotAngle[k])
x_store = append(x_store,x_val)
y_store = append(y_store,y_val)
ordered_x = concatenate((x_store,-x_store[::-1],-x_store,x_store[::-1]))
ordered_y = concatenate((y_store,y_store[::-1],-y_store,-y_store[::-1]))
straightened_cutout[j] = shapely_cutout(ordered_x,ordered_y)
minimizer_kwargs = {"args": (centred_cutout[j],straightened_cutout[j],), "method": 'Powell'}
x0 = np.array([-90])
# output = basinhopping(visualRotation,x0,
# niter=3,minimizer_kwargs=minimizer_kwargs)
# visual_angle[j] = output.x
visual_angle[j] = -45
straightened_cutout[j] = af.rotate(straightened_cutout[j], visual_angle[j])
In [210]:
# j = find(cutoutRef == 38)[0]
# i = cutoutRef[j]
# plot(centred_cutout[j].exterior.xy[0],centred_cutout[j].exterior.xy[1])
# plot(straightened_cutout[j].exterior.xy[0],straightened_cutout[j].exterior.xy[1],'r--')
# axis('equal')
# grid(True)
In [211]:
# zoneWeight = 3000*widthFit.derivative()(zoneWidth)**2
# plot(zoneWidth,zoneWeight)
In [212]:
# numZones = 1000
# zoneWidth = linspace(minVal,maxVal,numZones)
# zoneWeight = zoneWeight = 3000*widthFit.derivative()(zoneWidth)**2
# circle = [0,]*(numZones)
# circle[0] = sh.Point(0,0).buffer(minVal/2)
# zone = [0,]*numZones
# zone[0] = circle[0]
# for i in range(1,numZones-1):
# circle[i] = sh.Point(0,0).buffer(zoneWidth[i]/2)
# zone[i] = circle[i].difference(circle[i-1])
# zone[-1] = sh.Point(0,0).buffer(20).difference(circle[-2])
# zoneWidth = linspace(minVal,maxVal,numZones)
# zoneWeight = 70*widthFit.derivative()(zoneWidth)
# zoneWeight[zoneWeight < 0] = 0
# circle = [0,]*(numZones)
# circle[0] = sh.Point(0,0).buffer(minVal/2)
# zone = [0,]*numZones
# zone[0] = circle[0]
# for i in range(1,numZones-1):
# circle[i] = sh.Point(0,0).buffer(zoneWidth[i]/2)
# zone[i] = circle[i].difference(circle[i-1])
# zone[-1] = sh.Point(0,0).buffer(20).difference(circle[-2])
def weightedAreaCalc(inputShape):
weightedArea = 0.
for i in range(numZones):
weightedArea += inputShape.intersection(zone[i]).area * zoneWeight[i]
return weightedArea
def weightedAreaDifference(cutout,ellipse,a,b):
cutout = af.translate(cutout, xoff=-a, yoff=-b)
ellipse = af.translate(ellipse, xoff=-a, yoff=-b)
return (weightedAreaCalc(ellipse.difference(cutout)) +
weightedAreaCalc(cutout.difference(ellipse)))
# def weightedAreaDifference(cutout,ellipse,a,b):
# cutout = af.translate(cutout, xoff=-a, yoff=-b)
# ellipse = af.translate(ellipse, xoff=-a, yoff=-b)
# return ellipse.difference(cutout).area + cutout.difference(ellipse).area
def ellipse_in_cutout(minimiserInput,cutout,weight):
ellipse = create_ellipse(minimiserInput)
a = minimiserInput[0]
b = minimiserInput[1]
toMinimise = 100*weightedAreaDifference(cutout,ellipse,a,b)
return toMinimise
class MyTakeStep(object):
def __init__(self, stepsize=0.5):
self.stepsize = stepsize
def __call__(self, x):
s = self.stepsize
x[0] += np.random.uniform(-3*s, 3*s)
x[1] += np.random.uniform(-3*s, 3*s)
x[2] += np.random.uniform(-6*s, 6*s)
x[3] += np.random.uniform(-8*s, 8*s)
x[4] += np.random.uniform(-90*s, 90*s)
return x
def print_fun(x, f, accepted):
global successCount_basinhopping
global minVals_basinhopping
global nSuccess_basinhopping
a = x[0]
b = x[1]
w = abs(x[2])
l = abs(x[3])
if (l > w):
width = w
length = l
theta = np.mod(x[4],180)
else:
width = l
length = w
theta = np.mod(x[4]+90,180)
f = np.around(f,4)
if accepted:
if np.all(np.isnan(minVals_basinhopping)):
minVals_basinhopping[0] = f
successCount_basinhopping = 1
print(("first local minima of %.4f at:\n"+
"[a,b] = [%.4f, %.4f],"+
" [width, length] = [%.4f, %.4f],"+
" theta = %.4f\n")
% (f,a,b,width,length,theta))
elif (np.nanmin(minVals_basinhopping) == f):
minVals_basinhopping[successCount_basinhopping] = f
successCount_basinhopping += 1
print(("agreeing local minima of %.4f at:\n"+
"[a,b] = [%.4f, %.4f],"+
" [width, length] = [%.4f, %.4f],"+
" theta = %.4f\n")
% (f,a,b,width,length,theta))
elif (np.nanmin(minVals_basinhopping) > f):
minVals_basinhopping[0] = f
successCount_basinhopping = 1
print(("new local minima of %.4f at:\n"+
"[a,b] = [%.4f, %.4f],"+
" [width, length] = [%.4f, %.4f],"+
" theta = %.4f\n")
% (f,a,b,width,length,theta))
else:
print(("rejected local minima of %.4f at:\n"+
"[a,b] = [%.4f, %.4f],"+
" [width, length] = [%.4f, %.4f],"+
" theta = %.4f\n")
% (f,a,b,width,length,theta))
else:
print(("failed to find local minima at:\n"+
"[a,b] = [%.4f, %.4f],"+
" [width, length] = [%.4f, %.4f],"+
" theta = %.4f\n")
% (a,b,width,length,theta))
if successCount_basinhopping >= nSuccess_basinhopping:
return True
def straight_ellipse_calc(cutout,weight,userNSuccess,angle):
mytakestep = MyTakeStep()
global nSuccess_basinhopping
nSuccess_basinhopping = userNSuccess
global minVals_basinhopping
minVals_basinhopping = np.empty(nSuccess_basinhopping)
minVals_basinhopping[:] = np.nan
global successCount_basinhopping
successCount_basinhopping = 0
minimizer_kwargs = {"args": (cutout,weight), "method": 'L-BFGS-B',
"bounds": ((0,0),(0,0),
(minVal,maxVal),
(minVal,maxVal),
(angle,angle)) }
x0 = np.array([0,0,3,5,angle])
output = basinhopping(ellipse_in_cutout,x0,
niter=1000,minimizer_kwargs=minimizer_kwargs,
take_step=mytakestep, callback=print_fun)
return output
In [213]:
numZones
Out[213]:
In [217]:
ellipse_results = [0,]*len(cutoutRef)
# for j in range(len(cutoutRef)):
for j in [0]:
i = j
print("Cutout %.f" % (i))
ellipse_results[j] = straight_ellipse_calc(straightened_cutout[j],1,5,visual_angle[j])
print(" ")
In [218]:
ellipse_def = zeros([len(cutoutRef),5])
# for j in range(len(cutoutRef)):
for j in [0]:
ellipse_def[j,:] = array(ellipse_results[j].x)
In [219]:
# j = find(cutoutRef == 14)[0]
# i = cutoutRef[j]
ellipse = create_ellipse(ellipse_def[j,:])
scaled_fig_start(12,12)
# plot(centred_cutout[j].exterior.xy[0],centred_cutout[j].exterior.xy[1])
plot(straightened_cutout[j].exterior.xy[0],straightened_cutout[j].exterior.xy[1],'b')
plot(ellipse.exterior.xy[0],ellipse.exterior.xy[1],'g',linewidth=4,alpha=0.5)
text(-5,5,('Width: %.2f cm' % (ellipse_def[j,2])),fontsize=11)
text(-5,4.5,('Length: %.2f cm' % (ellipse_def[j,3])),fontsize=11)
scaled_fig_end(12,12)
In [220]:
# a = ellipse_def[:,0]
# b = ellipse_def[:,1]
# new_width = zeros([len(cutoutRef)])
# new_length = zeros([len(cutoutRef)])
# theta = zeros([len(cutoutRef)])
# w = abs(ellipse_def[:,2])
# l = abs(ellipse_def[:,3])
# ref = l > w
# new_width[ref] = w[ref]
# new_length[ref] = l[ref]
# theta[ref] = mod(ellipse_def[:,4],180)[ref]
# new_width[~ref] = l[~ref]
# new_length[~ref] = w[~ref]
# theta[~ref] = mod(ellipse_def[:,4]+90,180)[~ref]
In [221]:
# width[0:24] - new_width
In [222]:
# cutoutRef
In [223]:
# # Display table
# df = pd.DataFrame(transpose([cutoutRef, width[0:24], new_width, length[0:24], new_length, factor[0:24]]))
# df.columns = ['cutoutRef','width', 'new_width', 'length', 'new_length' , 'factor']
# df.to_csv(path_or_buf = "new_cutoutdata.csv", index=False)
# df
In [223]:
In [223]: